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Abstract 

In this paper, we develop a stable and fast numerical scheme for relativistic dissipative hydro- 
dynamics based on Israel-Stewart theory. Israel-Stewart theory is a stable and causal descrip- 
tion of dissipation in relativistic hydrodynamics although it includes relaxation process with the 
timescale for collision of constituent particles, which introduces stiff equations and makes prac- 
tical numerical calculation difficult. In our new scheme, we use Strang's splitting method, and 
use the piecewise exact solutions for solving the extremely short timescale problem. In addition, 
since we split the calculations into inviscid step and dissipative step, Riemann solver can be used 
for obtaining numerical flux for the inviscid step. The use of Riemann solver enables us to cap- 
ture shocks very accurately. Simple numerical examples are shown. The present scheme can be 
applied to various high energy phenomena of astrophysics and nuclear physics. 

Keywords: relativistic hydrodynamics; relativistic dissipation; relativity; kinetic theory; 
methods: numerical 



1. INTRODUCTION 

In recent years, various high energy astrophysical phenomena are extensively studied by us- 
ing the relativistic fluid approximation; for example, ultra relativistic jet |H Ql, GRB Jj, 0], 
Neutron star merger |]5l |6|], pulsar wind and accretion flows around massive compact 

objects ifjQlJ]. In addition, because of the recent finding of the strongly coupled quark-gluon 
plasma (sQGP) in the Relativistic Heavy-Ion Collider (RHIC), description by relativistic hydro- 
dynamics equations have been vigorously studied in the context of nuclear physics [11]. As is 
well known, relativistic fluid equations are highly nonlinear because of the Lorentz factor and 
enthalpy, and the high energy astrophysical phenomena are studied mainly by numerical simu- 
lation except for simplified cases. For this reason, various numerical formulations of relativistic 
fluid are investigated. However, most of the existing numerical schemes assume ideal fluid ap- 
proximation, and there are only few studies taking into account the dissipation lfl2i [l3i [Til [l5ll . 

Relativistic fluid equation can be obtained by tensor decomposition of the particle flux vector 
and energy-momentum tensor [16]. When one considers ideal fluid, those tensors are decom- 
posed by assuming homogeneity and isotropy in the comoving frame of each fluid element, and 
dissipation variables are defined as the deviation from the ideal part of those tensors. However, 
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relativistic fluid includes independent two characteristic directions, the particle flux vector and 
energy flux vector, and this results in uncertainty of the definition of the fluid 4-vector. For the 
direction of fluid 4-vector, Eckart IU7II adopted that of particle flux vector, and his decomposition 
is called Eckart formalism; Landau and Lifshitz 111 811 adopted that of energy flux, and their de- 
composition is called Landau-Lifshitz formalism. In addition to these well-known formalisms, 
various formalisms are proposed 

If one considers the dissipative fluid equations, dissipation variables have to be expressed by 
the fluid variables, density, pressure, and fluid velocity. Eckart, Landau, and Lifshitz presented 
as the expression of dissipation variables the relativistic extension of Navier-Stokes approxima- 
tion. However, it is well-known that the Navier-Stokes equation is parabolic partial differential 
equations. This means that the Navier-Stokes equation is acausal, and not appropriate for the 
relativistic equation. In addition, Hiscock and Lindblom 11221 12311 have shown that the relativis- 
tic Navier-Stokes equations includes unphysical exponentially growing modes, and unstable for 
small perturbation. For this problem, Israel and Stewart (12411 proposed new relativistic dissipa- 
tion fluid theory called Israel-Stewart theory, which takes into account the second-order deviation 
terms for entropy. This theory is hyperbolic equations, and has been shown to be causal and sta- 
ble. However, the equations of Israel-Stewart theory include 14 variables, and extremely complex 
in contrast to relativistic ideal fluid equations. In addition, the evolution equations of dissipation 
variables include parameters corresponding to relaxation timescales that are much shorter than 
hydrodynamical timescale, and the resultant stiff equations make it difficult to integrate numer- 
ically. For these reasons, applications of Israel-Stewart theory to physical problems are very 
limited. 

The first numerical scheme of Israel-Stewart theory is proposed by Molnar et al. (2010) 1 25] . 
They integrate fluid equation by using SHASTA as a shock capturing scheme, and integrate 
evolution equations of dissipation variables by using the ordinary explicit difference scheme. 
As explained above, evolution equations of dissipation variables are very stiff, and they have to 
use sufficiently high resolution to resolve both the macroscopic and relaxation timescale, which 
demands a exceedingly high numerical cost. 

We develop a new numerical scheme for relativistic dissipative hydrodynamics that can inte- 
grate accurately and efficiently. We split the fluid equations into inviscid part and dissipation part. 
The inviscid part corresponds to ideal relativistic fluid equations, and can be solved accurately 
by using a relativistic Riemann solver 11261 \27l l28l l29l l30ll3 ll l32[ l33l 13411 . The Riemann solver is 
a method to calculate numerical flux by using exact solution of the Riemann problems at the in- 
terfaces separating numerical grid cells, and can be used to describe the flows with strong shocks 
and sharp discontinuity stably and highly accurately. When the dissipation terms are small, the 
dynamics of fluid is dominated by the inviscid part, and we can obtain more accurate numerical 
results by means of the Riemann solver. As for the evolution equations of dissipation variables, 
we use the Piecewise Exact Solution method (PES) B35L l36l 13711 . PES is a numerical method 
for solving the stiff equation by using the formal solution. The use of formal solutions for those 
stiff equations eliminates the Courant condition for relaxation timescale. In this way, our new 
numerical scheme can describe relativistic dissipative hydrodynamic equations highly accurately 
and efficiently. Recently, similar methods were applied to solving resistive RMHD Il4ll5ll . 

This paper is organized as follows. In Sec. [2] we present the basic equations of relativistic 
hydrodynamics. We present explicit forms for relaxation timescale parameters near the equi- 
librium. In Sec. [3] detailed explanation of our new scheme is presented. In Sec. |4] analyses 
on stability and causality of relativistic dissipation hydrodynamics equations are presented by 
using simple scalar equation. In Sec. [5] some results of one-dimensional and multi-dimensional 
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simulations are presented. 



2. Basic Equations 

Throughout this paper, we use the units c — 1 , and Cartesian coordinates where the Minkowski 
metric tensor 77^ is given by rj^ = diag(-l, 1,1,1). Variables indicated by Greek letters take val- 
ues from to 3, and those indicated by Roman letters take values from 1 to 3. 

We define the convective time derivative D and the spatial gradient operator V ff as follows 

DA" 1 """ = tPA^-"", (1) 
Vo a«-*=^a«-«.. ( 2 ) 

where the tensor y y is a projection operator on the hyperplane normal to u M 

= jf v + „v (3) 



2.1. Ideal fluid 

The relativistic hydrodynamic equations can be obtained from the conservation of particle 
number, momentum, and energy. 



rpltV 

1 V 



0, 
0, 



(4) 
(5) 



where N M is the particle number density current and 7 Viv the energy-momentum tensor. When we 
consider the ideal fluid, they are given by 



N" = mf, 

T" v = pfufu v + prf, 



(6) 
(7) 



where n denotes the proper particle number density, p is the pressure, h = 1 + e + pip is the 
specific enthalpy, e is the specific internal energy, p = mn is the proper rest-mass density, and 
m is the rest-mass of the constituent particle, is the four-velocity of the fluid, satisfying the 
normalization condition: u^u^ = - 1 . 

When we use Cartesian coordinates, the evolution equations of a relativistic fluid are 

( D \ „ ( Dv j 



d_ 

dt 



m 
E 



dxi 



mV + pl'j 



m> 



0, 



(8) 



where v is the fluid three-velocity, D, m, E are the mass, momentum, and energy density relative 
to the laboratory frame, and /'■' is the unit tensor. In the laboratory frame, D, m, E are given by 



D = yp, 
m = phy 2 \, 
E = phy 2 - p, 



(9) 
(10) 
(11) 



where y = (1 - v 2 ) ^ 2 is the Lorentz factor. This is the most common form of perfect fluid 
equations for the numerical hydrodynamics. 
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2.2. Causal dissipative fluid 

When one considers relativistic dissipative fluid, the system is neither homogeneous nor 
isotropic, and the decomposition of the particle number density current N M and the energy mo- 
mentum tensor T^ v will change. In addition, since the dissipative fluid has two different charac- 
teristic direction iV and r 0/J , the definition of the four velocity u^ is generally not unique, that 
is, there is an uncertainty or freedom for our choice of the direction of u M . Since one decom- 
poses N M and T^ v with respect to u?, this means the form of relativistic dissipative fluid equation 
is not unique. In this paper, we consider only Eckart decomposition whose four velocity u 11 is 
parallel to W. For Landau-Lifshitz decomposition and other decompositions, see the following 
references |fl8[ [l9i |20[ |2lll . 

In the Eckart decomposition, the particle number density current N M and the energy-momentum 
tensor P' v are written as 



N" = nil", 

T" v = pfaiV + pif v + q"u v + q v u" + t» v 



(12) 
(13) 



where q M is the heat flux vector and t>" is the viscosity tensor. 

From Eqs. (HJl and (0, the evolution equations of relativistic dissipative fluid are given by 



d_ 

dt 



D ^ 

m' + q°u l + q'u° + t ' 
E + 2q°u° + r 00 



dxi 



mV + pP j + q'u j + q j u' + T ,j 
m j + q°u j + q j u° + T 0j 



= 0, 



(14) 



In contrast to the non-relativistic case, dissipation variables are differentiated with respect to time. 
For example, time derivative of the heat flux vector remains in the second line of Eq. (TFfl i even 
in the fluid rest frame. This is because the energy flux is identified with the momentum density. 
Then, if one uses the relativistic Navier-Stokes terms as the dissipative ones, the fluid equations 
becomes parabolic, and the characteristic velocity of this theory becomes infinity. This means 
that the dissipation variables evolve into equilibrium values within infinitely short time, and the 
time derivative of them diverges. Hiscock and Lindblom 11221 12311 proved that the relativistic 
Navier-Stokes theory adopting any definitions of four-velocity is unstable in the sense that small 
perturbation will diverge exponentially with time in any frame except for Landau-Lifshitz theory 
in its rest frame. In order to find stable and causal theory, Israel and Stewart developed the 
following second-order theory from relativistic Boltzmann equation 



DTI = —(Y1 NS -IT)- I n , 
D^=i« s -^")-/f. 



bet = -tf NS - <f) - 



(15) 
(16) 
(17) 



where and II are the shear viscosity and bulk viscosity defined as the traceless part and trace 
part of the viscosity tensor t /jv respectively, 



t^ v = ny v + if 



(18) 



T q , Tn, and t k are the relaxation times, and / is second-order terms, which are the product of 
dissipation variables and derivative of fluid variables. Note that they can be neglected in the 
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astrophysical application, since the gradient of fluid variables are not so steep. If one considers 
the application to the QGP, one can use the following abbreviations [25] 



/] 



n = in(v^ + Dln^), 

tf"u V + 7T' h 'u")Du A 

i^(W + mnf) 



where a> is the rotational part of u^ v defined as 



(19) 
(20) 

(21) 
(22) 

(23) 



When one adopts Landau-Lifshitz frame as is often the case with QGP applications, I q van- 
ishes B25I1 . 

The "second-order" means that the entropy current contains second-order terms in deviations 
from equilibrium 11381 12411 . The above equations take into account time derivative of dissipa- 
tion variables, that is, relaxation effect. Owing to these terms, the relativistic dissipative fluid 
equations become hyperbolic, and it allows the equations become stable and causal if one uses 
appropriate parameters. We discuss these appropriate parameters in Sec|4] 

The Navier-Stokes terms are given by 



<f N s = ~ K i" V + TuPu v,p) 



(24) 



(25) 



where k is the heat conduction coefficient, 77 is the shear viscosity coefficient, and £ is the bulk 
viscosity coefficient. The subscript NS means "Navier-Stokes" terms. If dissipation terms are 
small, we can rewrite Eq. (l24l > by substituting ideal equation of motion into u p u vp , and obtain 



(26) 



In this paper, we use Eq. ( 1261 1 as the Navier-Stokes heat flux vector. 

Heat flux vector has 4 components, and viscosity tensor 10 components. However, these are 
constrained by the following orthogonality conditions: 



q°u 



-q J Uj. 



(27) 
(28) 
(29) 



As a result, the number of physical degrees of freedom reduces to 3 and 6 respectively. 
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Eqs. ( fZSb and (f26t are given in the covariant form. When one considers the fiat Cartesian 
coordinate and one-dimensional problem, they reduce to 



TlidpUa- + do-Up) + U - -T] 



Tj{u x u'd t u' + (-1 + (u') 2 )d,u x + (1 + (u x ) 2 )d x u' + u'u x d x u x ) 



= - [77 {m-V<5 ; m' + (-1 + (u'fydiu 1 - + u^u x d x u l + u'^d^ 



(30) 



(31) 



T](d p u a + S^Mp) + I £ - - 77 ) gpc-dxu" 



= — [77 [u y u'd t u z + u z u'd t u y + u y u x d x u z + u z u x d x u y } 



^(dpKcr + do-Mp) + I £ - - 77 I gpa-d A U A 



^j]^u- L u'd,u x + u-^u'd/U 1 - + u ± u x d x u x + (1 + (m*) 2 )^" 1 } 



+ |f--J7|nV-< 



T ;vs _ 77 



= -T0 x u x -\e- -77 |(i + K) 2 )e, 



pn 



= -AC 



kV 1 - I-d t p) + (1 + (« x ) 2 ) (^7' - ^rd x p 



Ins = -J XM « 



ph 
T 

<9p - — dpp 
ph 



ph 



ph 

where 9 is the expansion of the fluid defined as 



u' \ d,T - —d t p ) + u x ( d x T - —j^xP 



= = d^uf + T%u a . 



(32) 



(33) 



(34) 



(35) 



(36) 



(37) 



In the above expressions, the viscosity tensor includes both shear viscosity and bulk viscosity. 
This is because this form is useful for the directional splitting explained in Sec. 13.41 
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2.3. Explicit forms of relaxation time 

Israel-Stewart theory includes some transport coefficients, that is, dissipation coefficients and 
relaxation times, and their values depend on distribution functions and the cross sections of 
collisions between constituent particles. Since we are interested in the fluid approximation, we 
consider only the cases close to the equilibrium |24 , 161. 

In this case, the explicit forms of relaxation time can be expressed as follows: 

rn = ^p T q = KTp u T K = 2rfit (38) 



where 



r/(T-l)=J3\l+5h//3-h 2 ), (39) 
a = (T - l)Q,**/TQp, o-i = -(T - l)Tp, (40) 

A.*" ft _(E^)' (4i, 



h 2 Q?p' ^ \ r / h P \ r-i 
^ = i + 6h/p fli = l+w-D/r 

Q = 3r - 5 + 3T/hfi, Q* = 5 - 3T + 3(10 - lY)h/p, (43) 
Q" = 5 - 3r + 3T 2 /(T - l)h 2 p 2 , (44) 



(3 — m/T, and h is the enthalpy. 

Asymptotic forms of relaxation time are: 



1. In the non-relativistic limit (J3 — > oo, F = 5/3) 



2^ 2 2^ tj 

Tn = , T Cj = — — , t„ = - (45) 

5 p 5p p 

2. In the ultra-relativistic limit (J3 -> 0, r = 4/3) 

m 5kT 3/? 

T n = t 9 = — , r ff = — (46) 

/3 4 p 4/? 2p 

The explicit forms of dissipation coefficients k, r\, £ depend on the particle interaction, and 
one can use appropriate dissipation coefficients for each problems. 



3. Numerical scheme 

In this section, we start description of our scheme for one-dimensional case. Multi-dimensional 
case is shown in Sec l3.4l 



3.1. Strang splitting method 

The relativistic dissipative fluid equation is the hyperbolic -relaxation one, and this has dif- 
ferent difficulty for the inviscid and dissipation part respectively. The difficulty of inviscid part 
of fluid equation results from the non-linearity of the fluid equation, and this exists even in non- 
relativistic ideal fluid equation; the difficulty of dissipation part is that the evolution equations of 
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dissipation Eqs. ( fT5l ). ( fT6b , and (T% are stiff equations. When one solves stiff equation by us- 
ing explicit differentiation, At must be shorter than relaxation timescale parameters for stability. 
However, relaxation timescale parameters rn, TV, t ? are generally much shorter than dynamical 
timescale of fluid, and it needs a heavy computational cost. 

In order to address these problems separately, we apply Strang splitting method and split the 
relativistic dissipative fluid equation as follows 



d_ 

at 



D 

m' 
E 



_d_ 

dxJ 



Dv J 
m'v j + pl'i 
m J 



(47) 



a_ 

at 



D 

tri + q°u' + q'u° + t Qi 
E + 2q°u° + r 00 



d_ 

dxi 



q'u J + q J u' + r J 
q°ui + q'u + T°J 



= 0, 



(48) 



First, the inviscid part Eq. (l47t can be solved accurately by using the Riemann solver. The 
Riemann solver is a method that calculates numerical flux by using exact or approximate solution 
of the Riemann problem at the cell boundary, and it is known that this method is stable and 
accurate. 

Next, we consider the dissipation part Eq. (l48l . The second terms of Eq. (l48l are the 
remaining part of the flux of Eq. (TT4l . and includes dissipation variables and four velocity u* 1 . 
For the second-order accuracy in time, one must use states evolved half time-step At/2. The first 
terms of Eq. (l48l are conserved variables, and they includes dissipation variables unlike non- 
relativistic case. In order to calculate stably, one has to substitute for these dissipation variables 
not the Navier-Stokes terms, but evolved ones by Eqs. (fl5t . (TTol i. ( TP7I ). Then, the first terms 
of Eq. d48l includes four velocity u 1 *. For this reason, one has to calculate them before the 
calculation of inviscid part, and save them until the dissipation part. 

In summary, we split the conserved variables U of Eq. ( I4B1 as follows 



U - U ideal + Udis 



(49) 



where 





' D * 




{ 




U ideal — 


m' 


s Udissip — 


q ^ + q'u + t 0/ 


(50) 




, E J 


2q°u° + t 00 





First, we calculate U = V ' ideal + Udissip- Then, we evolve U ideal by using Riemann solver. Next, 
we calculate IJ ideal + Udissip as the initial value of U of dissipation step, and integrate Eq. (l48l 
over the full time-step. U ideal is the ideal part of conserved variable evolved by using Riemann 
solver. By using the conservation-law form for Eqs. (l4Tt and d48b . this method satisfies the 
conservation law of mass, momentum, and energy within machine round-off error. 

In addition, we adopt q x , q y , q z , t 0x , t 0v , t 0z , r yz , r ZA , r"' for the primitive variables of the dis- 
sipation. This selection is based on the equality of the spatial direction, and is useful for the 
directional split. The other variables can be calculated by the orthogonality conditions Eqs. ( l27l i. 
d28l . and d29l . Since the dissipation variables are necessary for the primitive recovery proce- 
dure, we define the variables at the cell-center. For the calculation of the numerical flux at the 
cell-boundary, we use the average of the dissipation variables, and evolve timestep At/2 at the 
cell-boundary. 
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3.2. Stiff equation 

When one describes relativistic dissipative fluid by using the Israel -Stewart theory, one has 
to take into account the evolution of dissipation variables. However, those equations are stiff- 
equations, that is, the form of those equations are given by 

t U - ^ (51) 

' relax 

where T re i ax is the relaxation time, and this is the characteristic timescale of the evolution of U. If 
this timescale t is much shorter than the fluid timescale T/y„y, Eq. ( BTT l is called "stiff-equation", 
and the stability of an explicit scheme is achieved only with a timestep size At < T re i ax <k T/;„,y. 
In general, this constraint is more restrictive than the Courant-Friedrichs-Lewy (CFL) condition 
At < Ax/c c i, aracl , where c c i iarncr is the characteristic velocity of fluid, for example, sound velocity, 
Alfven velocity, and so on. This increases the computational cost exceedingly, which hinders 
the use of simple explicit scheme. For avoiding this timestep restriction, we apply the Strang- 
splitting technique. 

First, rewriting Eqs. (fT5t , (Q2), (TTTt in the coordinate dependent forms, they reduce to 

ijt +vJ l) n -i (nNs - n) -^ (52) 



rlz + ^-k^r^-^)-/,, (53) 



at ox> ) r n 
A^ +Vj ^¥ = T^-^-h. (54) 



Then, we split the above equations as follows 



and 



|-n = — (n NS - n), (58) 

at jT n 

|^v = J_ (7r Mv_^v ) (59) 

dt jT n iVi 

|^=— (<f NS -<f)- (60) 

dt JT q 



Eqs. ( T55l ). d56*b . and d57l i are the first-order advection equation with source terms, and one can 
solve them by upwind scheme. Eqs. d58l ). d59b . and d60b are stiff equation that requires special 
care. In our new method, we solve the stiff equations by using the PES method |35, 3a 371. 
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Since we use Strang-splitting technique, one can obtain the formal solutions of Eqs. (|58T >. (l59l >. 
and ( f60b as follows: 



n = (n -nA,j)exp 

^ = ( 7r V - 7 45) eX P 



f-fo 



t-t 



t-t 



(61) 
(62) 
(63) 



where subscript means the initial value. Since these are the formal solution, numerical calcu- 
lation of these terms remains stable irrespective of the timestep. In this way, the time-step of 
our scheme is not restricted by the stiff equations Eqs. (l58l ). d59l ), and d60b those correspond to 
relaxation. In this way, we solve the stiff equation by using piecewise exact solution Eqs. doTt - 
d63l . Thus, this procedure is called PES method. Note that the terms /; / in Eqs. dT9b and (1221 
include dissipation variables, which might complicate the actual exact solution. In this paper, we 
recommend to assume that the dissipation variables in these terms are constant, and simply add 
them to the Eqs. doTt - d63l . If this procedure results in bad approximation, it means that one 
should not use fluid approximation but the kinetic equation. This is because these terms should 
be small compare to the other dissipation terms when fluid approximation is justified. 

The Navier-Stokes terms Eqs. (l30l - (l36l include not only spatial derivatives but also time 
derivatives unlike non-relativistic case. We calculate time derivatives by using the following 
form of first-order explicit finite differentiation 



d,U" = 



u ideal 



jjn 



Ar 



(64) 



where U^}, is the variable evolved by the inviscid step, and U" is the initial value of the n-th 
step. For this reason, one has to save initial fluid variables until dissipation step. In addition, 
we approximate the spatial derivatives of the Navier-Stokes terms with the centered finite dif- 
ferences. This is because the physical meanings of the dissipation variables are the diffusion. 
For other part of the spatial derivatives, we use the MUSCL scheme by Van Leer M4011 for the 
second-order accuracy in space. 



3.3. Primitive recovery 

When one uses conservation-law form, updated variables are not primitive variables but con- 
served variables, and one has to obtain the primitive variables from conserved variables. In the 
case of relativistic fluid, one has to solve a non-linear algebraic equation for this primitive recov- 
ery even in the ideal fluid case. In the case of relativistic dissipative fluid, the primitive recovery 
is more complicated, since the conserved variables include dissipation terms. In this section, we 
explain a method that makes the primitive recovery somewhat simple. 

The dissipation variables in the conserved ones are obtained by using relaxation equations 
Eqs. ( fT5l ). ( TToT ). and (TTTb . Then, a major cause of the complexity of the primitive recovery is four 
velocity if multiplied by dissipative variables. If the effect of the dissipation is small enough, 
one can expect that the change of four velocity vf during dissipation step is smaller than during 
inviscid step. For this reason, we use four velocity if obtained after inviscid step as the initial 
guess, and this enable us to calculate the dissipation part of conserved variables Udissip- We can 
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obtain ideal part Ujdeai by subtracting Udissip from conserved variables U, and we can obtain 
primitive variables by using the primitive recovery for ideal fluid part. Then, we replace the 
initial guess for u 11 in Udissip with the obtained four velocity, and carry out the primitive recovery 
again. We repeat this procedure until the primitive variables converge. In this way, one can 
obtain primitive variables consistently. 
In summary, 

1. Calculate Udissip using vP obtained by the inviscid step as an initial guess, and evolve 
dissipative variables. 

2. Calculate U ideal = U — Udissip, and obtain primitive variables by using ordinary primitive 
recovery of relativistic ideal fluid. For the stability of the numerical calculation, we have 
imposed the stability condition \U\ > c\Udissip\, where c is some constant number smaller 
than unity. 

3. Replace the initial guess for u u in Udissip with the obtained four velocity, and carry out the 
primitive recovery again. 

4. Repeat this procedure until the primitive variables converge. 

From test calculations, this approach seems to work well even for large dissipation coeffi- 
cients and discontinuous profiles for physical variables, and converges within less than 5 itera- 
tions. 

3.4. multi-dimensional case 

We have explained one-dimensional scheme so far. For multidimensional calculation, we 
can apply the directional splitting method ll39ll where one applies one-dimensional operator in 
each spatial direction successively. To achieve second-order accuracy in time, one has to apply 
one-dimensional operator in the following order for two-dimensional case 

U" +l = l}J%l}J 2 U n , (65) 

and for three-dimensional case 

u n+x = Ly^^Ly^y^y^y^y 3 ^ 6 ^/ 3 ^/ 6 ^/ 3 ^^ 6 ^. (66) 

The At can be determined by Courant-Friederichs-Lewy (CFL) condition presented in the next 
section. 

4. Causality and Stability 

The Israel-Stewart theory is known as the stable and causal relativistic dissipative fluid the- 
ory. Strictly speaking, one has to use the appropriate parameters for the stability and causality. 
However, in order to guarantee the stability and causality, the values of the various dissipation 
coefficients should be limited to certain ranges. In this section, we discuss these parameters. 
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4.1. Stability of the telegrapher equation 

The essential structure of the Israel-Stewart theory can be analyzed by the following simple 
equations 

3,g + V-F = 0, (67) 
d t F=--(F + T,VQ). (68) 

T 

Eliminating F from the above equations, one obtains the following telegrapher equation 

dfQ + -d t Q - -AQ = 0. (69) 

T T 



Then, the characteristic velocity of Eq. ( 1691 is given by 



(70) 

For the causality, the characteristic velocity must be lower than the maximum physical velocity. 
Thus, the relaxation timescale parameter r has the following physical lower limit 

r > T min = -=— , (71) 

^max 

where v max is the maximum velocity of the physical system, that is smaller than the speed of 
light. 

Next, we consider the stability condition. Israel and Stewart [24] indicates that the Israel- 
Stewart theory is stable and causal theory if one considers the Boltzmann gas, that is, near the 
equilibrium distribution. When one uses a simple explicit finite difference scheme, the CFL 
condition is that the timestep size At is less than the relaxation time r. At < r. In the following, 
we consider the stability condition of PES method for Eqs. (|67| | and d68l by using the von 
Neumann's stability analysis and numerical simulation. 

In this case, the basic equations can be written as follows. 

d,Q + d x F = 0, (72) 
F = -rid x Q + (Fo + T,d x Q)e-T. (73) 

Since we use PES, the relaxation equation of F is replaced by the formal solution. Then, we 
differentiate Eqs. (l72l and (1731 . In the case of the second-order accuracy in time, 

O n+l — O" f" +l / 2 _ pn+l/2 

^7 7+1/2 /-1/2 „ 

= 0, (74) 



Af Ajc 

„ +l _ Uj+V2 Uj-w 
t j - -T] : + 



n n+l/2 _ n n+l/2 

' 1 Ax 



Ax 

We substitute into the above equations the following forms of Q" and F" 



(75) 



Q"j = R"e ijB , F 1 ] = G"e ij8 , 
12 



(76) 



where R n and G n are the n-th power of the constants R and G. Then, they reduce to 



At ft 
R n+l -R"+2 — /sin-G" +1/2 = 0, 
Ax 2 

G «+i = ——/sin —R n+l l 2 + G n e~^ + ^isinV'/VT. 
Ax 2 Ax 2 

From these equations, we can obtain the following forms of equation 



R" 



_Af At _Af 9 6 



R" + e--R"- x = 0. 



By solving this equation, we obtain 



R=- 

2 



1 +, 



At _i J 

-4—77(1 -e Osin z - 

Ax L 2 



*\4 



(77) 
(78) 

(79) 



, _ 41 . A? . -, At 

l+e t -4— 77 sin 2 -(1 -e ~) 
Ax z 2 



(80) 

For the stability, it is necessary to satisfy the criterion |/?| < 1 . In order to obtain the stability 
restriction, we substitute \R\ - 1 for Eq. (ISO! , and the equation reduces to 



2rj sin' 



— if At » t 

Af <sc r 



(81) 



( -sin f V-7 " 

The exact relation is given in Fig. Q] Eq. ( l8TT l and Fig. [T]show that the stability restriction of the 



100 



0.1 1 10 

A X / V I T) 

Figure 1: Numerical solution of Eq. {80} when |R| = 1. This figure shows that Af/r is proportional to Aa" when At/r <sc 1, 
and Ax 2 when Af/r » 1. 

upper limit for Af for our numerical scheme with PES is the same dependence on Ax as that of 
the parabolic equation At < Ax 2 /i] when Af > r, and it becomes the same dependence on Ax as 
that of the hyperbolic equation Af < Ax Jrjrj when Af < r. Q 



1 In addition, when one adopts first-order accuracy in time, the stability restriction becomes At < Ax 2 1 At] when 
At > T. 
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Next, we solve the telegrapher equation Eq. (l69l > numerically by using PES method for the 
relaxation equation Eq. (l68l . and study the dependence on the parameters. We take the number 
of grid points = 400, and use the CFL number 0.4. For the stability restriction, we use Eq. 
(ISTT i. The initial left and right states are given by 

Q L =3.0 for x < 1, (82) 

e R =1.0 for x>\. (83) 

We set the initial value of F as F = 0. 

First, we set the relaxation timescale t - 2.0 x 10~ 3 , and the dissipation coefficient r\ = 
1.0 x 10~ 3 . Fig. |2]is the numerical result of PES method, and it reproduces the diffusion of the 
Q very well. 

3.5 
3 
2.5 
O 2 
1.5 



0.5 

-1 -0.5 0.5 1 

x 

Figure 2: Numerical solution of telegrapher equation Eq. i69\ . A snapshot at t = 1.0 is shown. Diffusion smooths the 
initial discontinuity . For this test problem, a cell size Ax = 0.005 is used, and we set dt = 5.0 X 10~ 3 , r = 2.0 X 10~ 3 , 
and;; = 1.0 X 10~ 3 . 

Next, we consider the case of violating the stability restriction Eq. ( IBTb . In this case, we set 
the relaxation timescale t = 2.0 x 10~ 3 , the dissipation coefficient rj - 1.0 x 10~ 3 , and set the 
timestep At = 1.1Ax 2 /2tj. Fig. [3]is the numerical result, and it shows that the solution diverges 
when one violates the stability restriction obtained by the von Neumann's stability analysis. 

4.2. Stability conditions of new numerical scheme for Israel-Stewart theory 

The results of previous section will be able to be applied to the new numerical scheme for 
the Israel-Stewart theory, since the structure of the equations of the Israel-Stewart theory are 
telegrapher equations. For this reason, we impose the following stability conditions as the CFL 
condition 

I min \Ax/c s , ap/iAx 2 /max{/c, n, Z"}) if At > r, 
At=C a \ H (84) 

I min Ax/c.5, Ax yr/m if At < t 

where c s is the sound velocity, a is 1/2 when second-order accuracy in time and 1 /4 when first- 
order accuracy in time, and < C a < 1 is the Courant number. The above stability conditions 
are just provisional ones, since the dissipative RHD equations are highly non-linear equations, 
and it is difficult to obtain the exact conditions. Practically, the above conditions work well in 
our calculation. 
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Figure 3: Numerical solution of telegrapher equation Eq. d69t in the case of violating the stability condition (At > A(,„ m ). 
A snapshot at t = 0.25 is shown. The unstable modes grow quickly. For this test problem, a cell size A,t = 0.005 is used, 
and we set* = 1.375 X 10~ 2 , r = 2.0 X 10~ 3 ,and;7= 1.0 X 10~ 3 . 

Next, when one solves the Israel-Stewart theory by using conservation-law form, one has 
to recover primitive variables from conservative ones. Then, if one considers large dissipation 
coefficients, the dissipation part of conservative variables Udissip can be in general equal or larger 
than U ideal- This makes the primitive recovery unstable, since the error of dissipation variables 
affects considerably conserved variables U similar to the low f3 case of MHD equation. This may 
happen when the previous condition Ax < -\Jrrj is violated. This is because yjrri is equivalent 
to the mean free path in ideal gas, and Ax < s/rrj means that one resolves length scale shorter 
than the mean free path. Since the Israel-Stewart theory is approximation of the Boltzmann 
equation, the approximation becomes bad in this region. For this reason, one cannot use Israel- 
Stewart theory in such parameters. In contrast, if one considers effective theory, we recommend 
to impose the following restrictions on the dissipation coefficients 

q' < phy cs c s , (85) 
t ! U phy 2 cs c 2 s , (86) 

where c s is the sound velocity, and y cs = 1 / yl - cf. We recommend to impose this restriction 
during the evolution step of dissipation variables, the calculation of the numerical flux, and the 
primitive recovery procedure. Similar conditions can be found in the previous studies ll25ll . 

In addition, if one adopts t 0x , t°- v , t Q: , T yz , t zx , T xy for the primitive variables, one needs t xx 
for the numerical flux. In this case, if one calculates r xx by the orthogonality condition Eq. 
(1271 1. the numerical calculation sometimes becomes unstable, since the calculation includes the 
division by the velocity. In order to prevent this numerical divergence, we use the Navier-Stokes 
term for t xx . This approximation becomes bad when the timestep At is close to the relaxation 
time. However, if one considers the region where the fluid approximation is not so bad valid, the 
relaxation time is small and the above approximation gives accurate value for t xx . 

5. Test Calculations 

In this section, results of several one-dimensional and multi-dimensional test simulations are 
presented. In the following test problems, we use CFL number 0.4, and we consider an ideal 
equation of state h = 1 + r/(T - l)p/p. 
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5.1. ID test 



5.1.1. shearflow 

For the stability and causality, the Israel-Stewart theory adds 9 variables to the 5 fluid vari- 
ables, and it is very difficult to obtain exact solutions. However, relaxation of tangential velocity 
by the shear viscosity can be expected that the solution approaches to the exact solution of the 
non-relativistic case if one uses sufficiently low velocity. Thus, we consider the following initial 
condition 

(p L ,p L y L ) =(1-0,1.0,-0.1) for x<0.0, (87) 

(p R ,p s y R ) =(1.0,1.0,0.1) for x>0.0. (88) 

All the other fields are set to 0. 

In this case, the relativistic Euler equation is given by 

p/zd,[yV ] + dy y = 0. (89) 

Since we consider sufficiently low velocity, the Lorentz factor y is nearly unity. The viscous 
tensor r /iv reduces to 

r Ky = -jjflV. (90) 
Using Eqs. ( 1891 and j90l , the following exact solution can be obtained 




-0.15 1 1 1 1 1 

-1 -0.5 0.5 1 

x 



Figure 4: Numerical solution of the relaxation of a shear flow. Crosses denote a snapshot at t = 4.0. The analytical 
solution Eq. 19 It is shown as a solid curve. Our scheme for Israel-Stewart theory reproduces the analytical solution very 
well. For this test problem, a cell size Ax = 0.01 is used, and we set T = 4/3, and rj = 0.01. 

size A* = 0.01 is used, and the viscosity coefficient rj = 0.01. Fig. |4]shows that the numerical 
solution reproduces the exact solution very well. The convergence test shows that it is consistent 
with the second-order accuracy. 
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5.1.2. shock tube test 

In this section, we consider the shock tube problem. When one describes hydrodynam- 
ics including shock waves by using "ideal" fluid solver, the thickness of the discontinuities are 
determined by the numerical dissipation. In reality, the thickness of the shock wave front is 
determined by dissipation coefficients. In the following, we show that the thickness of the dis- 
continuities depends on the dissipation coefficients of our code. The exact Riemann Solver for 
ideal fluid part is based on the solution by Marti and Miiller [32]. 

We prescribe the initial left and right states 

(p L ,p L ,v yL ) =(10.0,10.0,0.2) for x < 0.5, (93) 

(p R ,p R ,v yR ) =(1.0, 1.0,-0.2) for x>0.5, (94) 

with an ideal equation of state with F = 5/3, and all the other fields are set to 0. Integration is 
carried until t = 0.4, and a cell size Ax = 0.0025 is used. The dissipation coefficients thermal 
conductivity k, shear viscosity 77, and bulk viscosity £ are assumed constant, and set 10~ 15 unless 
stated otherwise. 

First, Fig. [5] is the numerical results of the tangential velocity v y that changes the shear 
viscosity 77 = 0.01,0.05, 10~ 15 compared to the RHD exact solution. It shows that both contact 
discontinuity and shock are diffused, and the width of the discontinuity of 77 — 0.05 is about 5 
times greater than 77 = 0.01. 

0.25 

0.2 
0.15 

0.1 
0.05 


-0.05 

-0.1 
-0.15 

-0.2 

0.2 0.4 0.6 0.8 1 

x 

Figure 5: The numerical solution for the relativistic shock tube problem. Tangential velocity profile at t = 0.4 are shown 
by crosses in the cases of three different shear viscosity coefficients:^ = 10 -15 , 0.01, 0.05. The result of ideal gas is also 
shown by a solid line, fn this test problem, the widths of discontinuities are determined by the physical viscosity. For 
this test problem, a cell size Ax = 0.0025 is used, and we set T = 5/3, and the CFL number 0.4. 

Next, Fig. [6] is the numerical results of the temperature T = pip for various values of 
the thermal conductivity k = 0.01,0.05, 10~ 15 compared to the RHD exact solution. It shows 
that the heat flows from high temperature region to low temperature region, and both contact 
discontinuity and shock are smoothed, and the width of the discontinuity of k — 0.05 is about 5 
times greater than k = 0.01. 

Fig. [7] is the numerical results of the velocity v x for various value of the bulk viscosity 
f = 0.01,0.04, 10~ 15 compared to the RHD exact solution. The bulk viscosity diffuses the fluid 
expansion 9 = lit, and in the one-dimensional case the expansion reduces to 9 — d x u x . Fig. [7] 
shows that the gradient of v x in x-direction is smoothed out. 
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Figure 6: The numerical solution for the relativistic shock tube problem. Temperature profiles at t = 0.4 are shown by 
crosses in the cases of three different thermal conduction coefficients:*: = 10 - ' 5 , 0.01, 0.05. All the profiles are similar 
near the shock wave, but the width of contact discontinuity is determined by the thermal conductivity. For this test 
problem, a cell size A.v = 0.0025 is used, and we set T = 5/3, and the CFL number 0.4. 



5.2. Two-dimensional Kelvin-Helmholtz Instability 

In this section, we present the numerical results of the two-dimensional Kelvin-Helmholtz 
instability (KH instability) for the multidimensional numerical test problem. Multidimensional 
extension can be achieved simply via directional splitting explained in Sec. 13.41 

KH instability is that of a tangential discontinuity between parallel flows. It is well-known 
from the linear perturbation analyses that this instability arises when the Lorentz factor of fluid 
is not so high, and the growth rate is proportional to the wave number. Second property means 
that the numerical simulation of KH instability does not converge in the ideal fluid case; if one 
increases the number of grid points, grid-size-scale perturbations grow most rapidly, and this 
prevents numerical convergence of KH instability in the ideal fluid case. However, if one con- 
siders dissipation, the numerical convergence is possible since the perturbations shorter than the 
characteristic scale of dissipation are smoothed out. 

The initial condition is prescribed as 

(p, p, v\ v v ) =(1, 0.3, 0.1, 0.0) for y > 0.0, (95) 

(p,p,v x ,v y ) =(2,0.3,-0.1,0.0) for y < 0.0. (96) 

To trigger the KH instability, we perturb the shear flow by the position of tangential discontinuity 

y tangential 

y,angential = 0.01 Sm(k x X), (97) 

k x = 2n (98) 

We use an equation of state with F = 5/3, and the non-relativistic limit of relaxation time Eq. 
(l45l l. The computational domain covers the region [-1, 1] x [-1, 1] with 1024 x 1024 grid points. 
The CFL number is 0.2, and the integration is carried out until t = 30, that is, about 1.5 fluid 
crossing time. We set periodic boundary condition for x-direction, and reflecting boundary con- 
dition for y-direction for simplicity. 

Fig. |8]is the result of ideal fluid case carried by using a relativistic Godunov scheme l34ll for 
the sake of comparison. This figure shows that the rolling up of the interface results from the 
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Figure 7: The numerical solution for the relativistic shock tube problem. The profile of longitudinal velocity V at 
t = 0.4 is shown for three different shear viscosity coefficients = 10 -15 , 0.01, 0.05. The widths of discontinuities and 
rarefaction fronts are determined by the physical viscosity. For this test problem, a cell size A.v = 0.0025 is used, and we 
set T = 5/3, and the CFL number 0.4. 

KH instability. Note that numerical grid-size-scale perturbations grow in addition to the initial 
perturbation of sin function mode. Figs. |9]are the result of our new code of Israel-Stewart theory. 
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Figure 8: Non-linear development of Kelvin-Helmholtz instability in two-dimensional simulations without viscosity. 
The density profile at ; = 30 is presented. Numerical integration has been performed with Af = 0.2Ax/c s . Rolling up 
of the interface (tangential discontinuity) is the non-linear consequence of the rapid growth of the initial perturbation of 
"sine" function, characteristics of the KH instability. However, numerical cell-size-scale perturbations also grow, and 
contaminate the result. 

In this simulation, we consider only shear viscosity, and use r\ = 0.005 and r\ = 0.001 as the 
dissipation coefficient. Figs. |9]show that numerical grid-size-scale perturbations are stabilized, 
and only the initial perturbation of sine mode forms vortices. Note that for avoiding the numerical 
grid-size-scale perturbation in inviscid case, one has to introduce a sufficiently large scale of 
gradient a to the profile of fluid variables d, v x as follows: 

e = e + Aetanh(v/a), (99) 



where Q = (Q y>0 + Q y<0 )/2 and AQ = \Q y>0 - Q y<0 \. 

19 



0.300O0000O00000O00D+O2 




1 



Figure 9: Non-linear development of Kelvin-Helmholtz instability in two-dimensional simulation with viscosity. The 
density profile at t = 30 is presented. Numerical integration has been performed with At = 0.2A/„„„, where Ar,„,„ is 
presented in Eq. ( 18 It . We set the shear viscosity 1/ = 5 X 10~ 3 for the upper fluid, and rj = 10T 3 for the lower fluid. 
Numerical cell-size-scale perturbations are stabilized by the viscosity. 

Fig. [10]is L\ norm errors of the density under different grid points carried out by the new 
dissipation code in the case of 77 = 0.005. The solutions are compared to the result of 1024 grid 
points. This figure shows that the numerical results converge to the result of 1024 grid points 
with the second-order accuracy because of the shear viscosity. 

In Fig. QT| we compare the time evolution of the amplitude of v y that is a square root of 
the sum total of to the analytical solution of the linear perturbation obtained by Turland and 
Scheuer 114 ill . The numerical integration is carried out by the new dissipation code in the case 
of 77 = 0.005,0.001, and inviscid case. During the initial phase, the evolution of v v reproduces 
the prediction of the linear theory well. After the linear phase, the numerical solutions grow 
non-linearly, and start to deviate from linear theory because of the effect of shear viscosity and 
numerical grid-size-scale perturbation. Note that when the shear viscosity 77 takes larger value, 
the growth of v y saturates faster, and decay in time. 
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Figure 10: L\ norm errors of the density of the 2-dimensional shear flow calculated by the new dissipation code in the 
case of r\ = 0.001. This figure shows that the numerical solutions converge to the result of 1024 grid points with the 
second-order accuracy. 

6. Conclusion 

In this paper, we have presented a new numerical scheme for relativistic dissipative hydro- 
dynamics, that is, Israel-Stewart theory. Israel-Stewart theory is a stable and causal relativistic 
dissipation theory. However, for the stability and causality, this theory includes relaxation equa- 
tions of dissipation variables. In general, relaxation timescales of dissipation variables are much 
shorter than characteristic timescale of hydrodynamics. This means that relaxation equations 
of dissipation variables are stiff equations, and this makes it difficult to integrate Israel-Stewart 
theory numerically. In our new scheme, we use Strang's splitting method, and obtain formal 
solution of the relaxation equation for solving this extremely short timescale problem. By us- 
ing the formal solution, the Courant condition of relaxation time disappear, since we do not use 
explicit finite difference scheme. In addition, since we split the calculation of inviscid step and 
dissipation step, Riemann solver can be used for obtaining numerical flux of inviscid part, and 
this enables us to obtain more accurate solution. 

In astrophysical application, it is very important to take into account dissipation terms, since 
dissipation terms transform kinetic energy of bulk fluid into thermal energy, which becomes 
observable as thermal radiation. In addition, if one considers the dynamics of accretion disk, a 
viscosity is often used for the phenomenological models of angular momentum transfer, and our 
new scheme can be used for those modeling. In recent years, the strongly coupled quark-gluon 
plasma (QGP) in the Relativistic Heavy-Ion Collider (RHIC) has been vigorously studied by 
using description by relativistic dissipative hydrodynamics equations in nuclear physics, and our 
scheme will be useful for such calculations. 
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Figure 1 1 : Evolution of the amplitude of v y as a function of time comparing to the analytical solution of the linear 
perturbation obtained by Turland and Scheuer 114 ill . The numerical integration is calculated by the new dissipation code 
in the case of rj = 0.005, 0.001, and inviscid case. This figure shows that the numerical solutions reproduce the analytical 
theory in the linear growth region. 
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